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ABSTRACT 


Emitter  localization  is  a  very  important  communications  tool  that  will  be 
extremely  valuable  to  a  multitude  of  different  military  as  well  as  civilian  applications. 
In  many  parts  of  the  world,  GSM  is  the  preferred  method  of  modulation  used  in 
mobile  phone  traffic.  This  thesis  addresses  the  time  difference  of  arrival  estimation 
applied  to  GSM  type  signals  using  wavelet-based  techniques.  Signals  are  generated 
using  the  Hewlett-Packard  Advanced  Design  System  software  and  processed  using 
algorithms  based  on  Matlab.  The  results  of  this  thesis  prove  that  we  can  improve  upon 
the  localization  of  a  GSM  emitter  through  the  use  of  wavelet-based  denoising 
techniques. 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  ability  to  obtain  localization  information  from  a  cellular  phone  transmission 
is  increasingly  becoming  a  much  sought-after  facility.  Thanks  to  some  improved 
techniques  in  cellular  communications,  this  ability  is  very  much  realizable.  As  a  matter 
of  fact,  recently  there  has  been  legislation  which  requires  the  development  of  this 
technology.  When  a  911  call  is  made  from  a  wire  line  phone,  the  dispatcher 
automatically  can  provide  the  location  of  where  the  call  was  made.  Currently,  there  is  no 
such  capability  for  calls  made  from  cellular  phones.  In  October  1999,  President  Clinton 
signed  a  law  which  makes  911  the  official  nationwide  emergency  number  for  both 
regular  and  cellular  phones  [Appendix  A].  This  legislation  followed  a  move  by  the  FCC 
which  requires,  by  the  end  of  2001,  that  cellular  911  calls  automatically  provide  the 
caller’s  location  [Appendix  A]. 

There  are  countless  ways  in  which  the  development  of  this  technology  will  be 
beneficial  to  many  different  public  and  private  enterprises.  If  this  technology  is  utilized 
to  its  full  potential,  it  will  completely  change  the  way  the  entire  wireless  industry  does 
business  as  well  provide  an  invaluable  resource  for  many  military  applications.  Listed 
here  are  just  a  few  of  the  benefits: 

■  Location  Sensitive  Billing  -  Network  operators  will  be  able  apply  charges 
depending  on  the  location  of  the  mobile.  In  this  way,  different  rates  can 
be  charged  for  calls  from  the  home  or  office.  This  would  allow  network 
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operators  without  copper-based  public  switched  telephone  network 
(PSTN)  to  offer  competitive  rates  for  calls  from  the  home  or  office.[l] 
Intelligent  Transport  Systems  -  Traffic  congestion  can  be  determined  by 
the  number  of  cellular  calls  coming  from  a  particular  stretch  of  highway. 
Vehicles  with  installed  maps  can  provide  directions  and  alternate  routes  to 
travelers  to  avoid  congestion  and  traffic  mishaps.  Also  companies  could 
track  the  location  of  all  the  vehicles  in  its  fleet  and  more  efficiently 
dispatch  vehicles  according  to  the  traffic  situation  and  distance  to  the  final 
destination. 

Enhanced  Network  Performance  -  By  monitoring  the  frequency  of  calls 
from  various  locations  over  time,  cellular  networks  can  be  better 
organized  to  provide  the  optimum  service  for  a  particular  coverage  area. 
This  includes  the  ability  to  accurately  monitor  the  movement  of  the  mobile 
telephone  to  enable  the  network  to  make  better  decisions  on  when  to 
initiate  call  handoff  to  a  different  cell. 

Electronic  Surveillance  -  Law  enforcement  authorities  can  use  position 
information  to  track  the  location  of  a  mobile  station  which  could  serve  a 
multitude  of  purposes.  [2] 

Automobile  Security  -  Location  information  can  be  gained  from  a 
strategically  placed  mobile  transmitter  in  the  event  of  a  car  theft.  This 
would  aid  tremendously  to  the  recovery  of  stolen  vehicles  and  greater  rate 
of  apprehension  of  car  thieves. 


■  Unauthorized  Cellular  Usage  Detection  -  Cellular  carriers  have  long  been 
plagued  by  the  unauthorized  access  to  personal  communications  systems  . 
These  illegal  activities  would  diminish  drastically  with  the  advent  of 
cellular  phone  localization. 

B.  OBJECTIVE 

There  are  several  methods  available  that  would  allow  us  to  implement  a  position 
location  capability  into  cellular  communications  systems.  For  reasons  that  will  become 
more  apparent  in  the  following  chapters,  the  Time  Difference  of  Arrival  (TDOA)  method 
appears  to  be  the  most  logical  choice.  In  this  thesis,  we  will  discuss  the  TDOA 
estimation  and  its  application  to  Global  System  for  Mobile  (GSM)  signals.  Specifically, 
we  will  investigate  several  methods  to  reduce  the  mean  squared  error  of  the  TDOA 
estimate,  thereby  allowing  a  position  localization  with  a  corresponding  lower  error 
probablity. 

GSM  was  chosen  to  provide  the  backbone  for  this  research  because  of  its 
worldwide  acceptance  as  the  Personal  Communication  Service  (PCS)  network  of  choice 
and  its  inherent  positioning  localization  advantages.  We  will  discuss  this  system  in  great 
detail  in  order  to  provide  an  understanding  of  its  advantages.  We  will  also  demonstrate 
our  ability  to  simulate  GSM  signals  using  the  Hewlett-Packard  Advanced  Design  System 
(HP-ADS)  program.  Finally,  these  generated  GSM  signals  will  be  evaluated  using 
Matlab-based  programs  to  provide  an  indication  of  how  to  improve  upon  the  localization 
of  a  GSM  emitter. 
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C.  RELATED  WORK 


The  initial  research  on  this  project  was  begun  in  1998  by  Unal  Aktas,  LTJG, 
Turkish  Navy.  His  work  focused  on  the  use  of  the  wavelet  transform  to  increase  the 
accuracy  of  the  TDOA  estimation.  He  investigated  several  different  denoising  techniques 
and  presented  a  comparison  of  the  techniques.  These  included  wavelet  denoising  based 
on  Donoho’s  method,  wavelet  denoising  using  the  Wo-So-Ching  threshold,  wavelet 
denoising  using  hyperbolic  shrinkage,  wavelet  denoising  using  median  filtering, 
modified  approximate  maximum-likelihood  delay  estimation,  denoising  based  on  the 
fourth  order  moment,  and  a  time-varying  approximate  maximum-likelihood  technique 
The  results  of  his  work  showed  that  the  fourth  order  moment,  and  the  modified 
approximate  maximum-likelihood  techniques  were  superior  to  the  others.  For  that 
reason,  this  thesis  only  makes  use  of  those  three  techniques.  Modified  versions  of  the 
Matlab  code  for  wavelet-based  denoising  developed  by  Aktas,  are  used  in  this  thesis. 
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II.  GSM  AND  THE  TIME  DIFFERENCE  OF  ARRIVAL 


A.  REVIEW  OF  GSM 

Global  System  for  Mobile  Communications  (GSM)  is  the  European  standard  for 
digital  mobile  telephones.  It  was  developed  in  1990  and  first  implemented  in  1991  to 
replace  first  generation  European  cellular  systems,  which  were  generally  incompatible 
with  each  other.  It  has  gained  worldwide  acceptance  as  the  first  universal  digital  cellular 
system  with  modem  network  features  extended  to  each  mobile  user,  and  is  a  strong 
contender  for  Personal  Communication  Service  (PCS)  above  1 800  MHz  throughout  the 
world.  [3] 


1.  GSM  Channel  Structure 

One  characteristic  of  GSM  that  makes  it  so  desirable  is  that  it  uses  both  Time 
Division  Multiple  Access  (TDMA)  and  Frequency  Division  Multiple  Access  (FDMA)  to 
transmit  and  receive  information.  TDMA  allows  for  multiple  conversations  Jo  take  place 
simultaneously  at  the  same  frequency  (channel)  by  rationing  different  time  slots  to  each 
user.  GSM  uses  two  25  MHz  bands  of  the  radio  frequency  spectrum,  one  for  uplink 
(mobile  station  to  base  station)  and  the  other  for  downlink  (base  station  to  mobile). 
FDMA  divides  each  of  these  blocks  into  125  channels,  each  with  a  bandwidth  of  200 
kHz.  TDMA  allows  for  8  different  time  slots  ,  as  shown  in  Figure  2.1  from  Ref  [1], 
with  each  time  slot  being  576.92  microseconds  in  duration.  The  grouping  of  these  8  time 
slots  is  called  a  frame.  This  combined  with  the  mechanics  of  FDMA  provides  for  a  total 
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of  1000  actual  speech  or  data  channels.  Currently,  there  is  a  half-rate  speech  coder  in  the 
process  of  development  that  will  double  the  number  of  channels  to  2000  [4], 


as  control  channels  or  logical  channels.  Each  logical  channel  has  a  specific  purpose 
ranging  from  transmitting  system  parameters,  to  coordinating  channel  usage  between 
base  station  and  mobile.  The  logical  channels  are  mapped  onto  dedicated  time  slots  of 
particular  frames.  This  is  executed  using  a  burst  structure,  the  most  commonly  used 
being  a  normal  burst.  Contained  in  the  middle  of  this  burst  is  a  26-bit  training  sequence, 
chosen  for  its  correlation  properties.  This  received  burst  is  correlated  with  a  local  copy  of 
the  training  sequence  to  allow  for  the  estimation  of  the  impulse  response  of  the  radio 
channel,  which  aids  in  the  demodulation  of  the  bits  in  the  burst.  Another  important 
characteristic  of  training  sequences  in  the  different  burst  structures  of  GSM  is  that  they 
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lend  themselves  to  the  measurement  of  time  which  can  be  used  to  aid  in  the  gathering  of 
time-based  positioning  information  [1]. 

2.  GMSK 

GSM  uses  0.3  Gaussian  Minimum  Shift  Keying  (GMSK)  as  its  method  of 
modulation.  This  facilitates  the  use  of  narrow  bandwidth  and  coherent  detection 
capability.  In  GMSK,  the  rectangular  pulses  are  passed  through  a  Gaussian  filter  prior  to 
their  passing  through  a  modulator.  Within  this  filter,  baseband  Gaussian  pulse  shaping 
smoothes  the  phase  trajectory  of  the  MSK  signal,  and  thus  stabilizes  the  instantaneous 
frequency  variations  over  time.  The  result  of  this  is  significantly  reduced  sidelobe  levels 
in  the  transmitted  spectrum.  The  value  0.3  represents  the  3  dB-bandwidth-bit  duration 
product  (BT)  or  normalized  pre-Gaussian  bandwidth,  which  corresponds  to  a  filter 
bandwidth  of  81 .25  kHz  for  an  aggregate  data  rate  of  270.8  Kbps.  As  shown  in  Figure 
2.2  from  Ref.  [4],  sidelobe  levels  fall  off  dramatically  with  decreasing  BT  product. 
However,  reducing  the  BT  product  increases  the  irreducible  error  rate  which  is  a  product 
of  intersymbol  interference  (ISI)  created  in  the  Gaussian  filtering  process.  GMSK 
sacrifices  the  irreducible  error  rate  for  extremely  good  spectral  efficiency  and  constant 
envelope  properties.  This  is  a  good  tradeoff  as  long  as  the  irreducible  error  rate  is  less 
than  that  produced  by  the  mobile  channel. 

All  base  stations  transmit  a  GSM  synchronization  burst  on  the  synchronization 
channel.  Contained  in  the  burst  is  another  useful  64-bit  training  sequence.  The 
autocorrelation  function  for  this  sequence  reveals  a  main  correlation  peak  with  a  width  of 
4  bits.  Compared  to  signals  used  in  other  time-based  positioning  systems,  this  is 
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relatively  wide.  The  GMSK  modulation  technique  is  responsible  for  this  increased  width. 
This  is  significant  from  a  time-based  positioning  standpoint  because  it  reduces  the  timing 
resolution  and  hence  the  maximum  achievable  position  accuracy. 
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3. 


Interference 


In  the  mobile  radio  path,  the  GSM  signal  will  experience  some  I  SI.  The  mobile 
radio  channel  is  generally  hostile  in  nature  and  can  affect  the  signal  to  such  an  extent  that 
the  system  performance  is  seriously  degraded  [5].  Inherent  in  any  communication  system 
is  the  phenomenon  known  as  multipath.  This  refers  to  the  situation  where  a  signal 
propagating  from  a  transmitter  to  a  receiver  can  travel  via  several  different  paths 
including  line  of  sight  and  one  or  more  reflected  paths.  The  reflected  paths  cause  delay, 
attenuation  in  the  amplitude,  and  phase  shifts  relative  to  the  direct  path.  In  some 
environments,  multipath  components  may  be  delayed  by  up  to  30  microseconds.  These 
delayed  components  cause  ISI  in  the  received  signal.  [6]  An  example  of  this  interference 
is  illustrated  in  Figure  2.3  from  Ref.  [1].  As  shown,  it  is  a  distorted  correlation  peak 
which  can  cause  errors  in  any  time  measurements.  The  degree  of  distortion  depends  on 
the  relative  amplitude,  delay,  and  phase  of  the  received  signals. 

The  conventional  GSM  receiver  contains  a  signal  demodulator  which  equalizes 
the  received  signal  to  improve  system  performance  in  the  presence  of  multipath  by 
combining  the  information  received  from  the  different  arrivals.  This  process  is  called 
adaptive  equalization.  The  equalizer  does  not  reject  multipath  signals,  but  combines 
them.  However,  for  positioning  reasons,  the  receiver  needs  to  be  able  to  determine  the 
first  arrival  corresponding  to  the  line  of  sight,  and  thus  requires  a  multipath  rejection 
algorithm  [1]. 
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Figure  2.3:  Comparison  of  Ideal  and  Multipath  Distorted  Correlation  Functions  for  GSM 
Synchronization  Burst  (64-Bit  Training  Sequence) 

4.  GSM  Coverage  Area  Layout 

The  area  that  a  GSM  network  covers  is  divided  into  a  number  of  cells,  each 
served  by  its  own  base  station  (BS),  also  known  as  a  base  transceiver  station.  The 
number  of  cells  served  by  each  BS,  called  a  cluster,  is  determined  by  the  amount  of 
system  capacity  (total  available  channels)  required  for  a  particular  area.  A  larger  cluster 
size  is  an  indication  of  a  larger  capacity  system.  One  tradeoff  associated  with  the  larger 
capacity  system  is  that  a  more  complex  system  is  required  to  control  the  increased 
demand  for  frequency  reuse.  Also,  more  cells  per  cluster  translates  into  smaller 
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individual  cell  size  and  thus,  more  co-channel  interference  for  which  to  account,  as  well 
as  a  greater  frequency  of  handoffs  between  cells.  The  typical  cluster  size  is  4,  7,  or  12. 

Each  cell  is  distinguished  by  a  unique  cell  identifier  and  is  allocated  one  or  more 
uplink/downlink  frequency  pairs.  Multiple  BS’s  may  be  grouped  together  under  the 
control  of  a  base  station  controller  (BSC).  In  turn,  several  BSC’s  are  usually  controlled 
by  a  mobile  service  center  (MSC),  which  handles  tasks  such  as  call  routing  and  serves  as 
the  interface  between  the  mobile  network  and  the  fixed  telephone  network. 

5.  GSM  Performance  Enhancing  Techniques 

In  addition  to  increasing  the  number  of  cells  per  cluster,  GSM  uses  a  number  of 
other  methods  to  increase  capacity.  These  are  important  because  they  provide  favorable 
implications  for  positioning.  One  such  technique  is  the  use  of  sectored  cells.  Sectoring 
effectively  breaks  the  cell  down  into  groups  of  frequencies  (channels)  that  are  only  used 
within  a  particular  region.  This  is  best  exemplified  as  slicing  a  pie  into  equal  segments. 
Not  only  does  this  technique  offer  additional  information  for  use  by  a  positioning  system, 
but  it  also  has  the  effect  of  reducing  co-channel  interference.  This  reduced  co-channel 
interference  is  due  to  the  fact  that  a  given  cell  receives  interference  and  transmits  with 
only  a  fraction  of  the  available  co-channel  cells. 

Another  performance  enhancing  technique  employed  by  GSM  is  slow  frequency 
hopping;  so  termed  because  the  hopping  rate  is  low  compared  with  the  symbol  or  bit  rate. 
The  mobile  radio  channel  is  a  frequency-selective  fading  channel,  which  means  the 
propagation  conditions  are  different  for  each  individual  radio  frequency.  Slow-frequency 
hopping  is  used  instead  of  fast-frequency  hopping  because,  in  order  to  perform  well,  a 
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frequency  synthesizer  must  be  able  to  change  its  frequency  and  settle  quietly  on  a  new 
one  within  a  fraction  of  one  time  slot  (576.92  microseconds).  Frequency  hopping 
provides  frequency  diversity  to  overcome  Rayleigh  fading  due  to  the  multipath 
propagation  environment.  Frequency  hopping  also  provides  interference  diversity.  The 
amount  of  interference  varies  with  any  given  channel  within  any  given  cell.  A  receiver 
set  to  a  channel  with  strong  interference  will  suffer  excessive  errors  over  long  strings  of 
bursts.  Frequency  hopping  prevents  the  receiver  from  spending  successive  bursts  on  the 
same  high-interference  channel.  These  beneficial  effects  of  frequency  hopping 
contribute  to  a  reduction  in  the  minimum  signal-to-noise  ratio  (SNR)  required  for  good 
communications. 

B.  POSITIONING  TECHNIQUES 

It  is  important  to  remember  that  in  order  to  comply  with  the  FCC  ruling,  any 
positioning  solution  must  work  with  existing  phones.  For  this  reason,  most  cellular 
geolocation  proposals  are  multilateral,  where  the  estimate  of  the  mobile’s  position  is 
formed  by  a  network  rather  than  by  the  mobile  itself.  Hence,  the  unilateral  proposal  of 
adding  a  global  positioning  system  (GPS)  capability  to  mobile  phones  is  not  a  universal 
solution  since  mobile  telephone  operators  would  have  to  replace  or  retrofit  every  mobile 
telephone.  Even  if  it  were  feasible  to  carry  out  this  retrofit,  there  are  other  factors  that 
would  make  this  a  futile  attempt.  First  of  all,  the  size  and  cost  of  mobile  transmitters 
would  be  considerably  increased  because  of  the  incorporated  GPS  receiver.  Moreover,  a 
clear  view  of  the  sky  is  necessary  to  receive  usable  GPS  data  from  at  least  three  satellites, 
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which  makes  the  system  useless  in  areas  surrounded  by  buildings,  mountains,  and  other 
obstructions — common  working  environments  for  cellular  communications.  [7] 

There  are  three  ways  that  position  localization  information  can  be  determined 
using  the  signaling  aspects  inherent  in  the  GSM  specification.  One  is  through  the  Angle 
of  Arrival  (AOA)  method,  which  uses  sector  information.  Second,  is  through  the  use  of 
propagation  time  using  timing  advances  (TA).  Third,  is  the  Time  Difference  of  Arrival 
(TDOA)  method,  which  offers  several  advantages  over  the  previous  methods.  In  addition 
to  providing  some  immunity  against  timing  errors  when  the  source  of  major  reflections  is 
near  the  mobile,  the  TDOA  method  is  less  expensive  to  put  in  place  than  the  AOA 
method.  Also,  the  AOA  and  TOA  methods  are  hampered  by  requirements  for  line  of 
sight  (LOS)  signal  components,  whereas  the  TDOA  method  may  work  accurately  without 
a  LOS  component  [8].  Each  method  will  be  briefly  described  before  taking  an  in-depth 
look  at  the  TDOA 

1.  Angle  of  Arrival  (AOA) 

If  two  or  more  BS  platforms  are  used,  the  location  of  the  desired  mobile  in  two 
dimensions  can  be  determined  from  the  intersection  of  two  or  more  lines  of  bearing 
(LOB),  with  each  LOB  being  formed  by  a  radial  from  a  base  station  to  the  mobile. 
Figure  2.4  illustrates  the  method  known  as  triangulation  using  three  different  platforms. 
Here,  the  LOB’s  from  the  mobile  to  two  adjacent  platforms  form  what  is  called  a 
measured  angle,  <}>.  After  two  measured  angles  are  calculated,  trigonometry  or  analytic 
geometry  may  then  be  used  to  deduce  the  location  of  the  MS. 
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Figure  2.4:  Triangulation 

Another  variation  of  the  AOA  method  relies  upon  the  sectoring  of  cells  for  all  of 
its  positioning  information.  Cell  sectoring  involves  the  use  of  highly  directive  antennas 
to  effectively  divide  the  cell  into  multiple  sections.  The  angle  of  arrival  can  be 
determined  at  the  base  station  by  electronically  steering  the  main  lobe  of  an  adaptive 
phased  array  antenna  in  the  direction  of  the  arriving  mobile  signal.  In  this  case,  a  single 
platform  may  be  sufficient  for  AOA  positioning,  although  typically,  two  closely  spaced 
antenna  arrays  are  used  to  dither  about  the  exact  direction  of  the  peak  incoming  energy  to 
provide  a  higher  resolution  measurement.  [9] 

AOA  measurements  might  provide  an  excellent  solution  for  wireless  transmitter 
localization  were  it  not  for  the  fact  that  this  method  requires  signals  coming  from  the 
mobile  to  the  base  station  be  from  the  line  of  sight  (LOS)  direction  only.  This  is  seldom 
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the  case  in  cellular  systems  given  that  the  operating  environment  is  often  heavily 
cluttered  with  rough  terrain,  tall  buildings  and  other  obstructions. 

2.  Time  of  Arrival  (TOA) 

Since  GSM  is  a  TDMA  system,  its  successful  operation  depends  on  the  ability  of 
all  signals  to  arrive  at  the  BS  at  the  appropriate  time.  And  since  the  signals  arriving  at  the 
BS’s  originate  from  different  distances,  the  time  at  which  the  signals  are  sent  must  be 
varied.  This  is  accomplished  by  having  each  BS  send  a  timing  advance  (TA)  to  each  MS 
connected  to  it.  The  TA  is  the  amount  by  which  the  mobile  station  (MS)  must  advance 
the  timing  of  its  transmission  to  ensure  that  it  arrives  in  the  correct  time  slot.  Obviously, 
the  TA  is  a  measure  of  the  propagation  time,  which  is  proportional  to  the  distance  from 
the  BS  to  the  MS.  Hence,  the  location  of  the  MS  can  be  constrained  to  a  circular  locus 
centered  on  the  BS.  If  the  MS  could  somehow  be  forced  to  hand  over  to  two  more  BS’s, 
it  would  be  possible  to  implement  a  positioning  scheme  based  on  the  intersection  of  the 
three  circles. 

One  problem  with  the  TOA  technique  is  that  artificially  forced  handoffs  can  be 
made  to  less  optimal  BS’s.  This  degrades  call  quality  as  well  as  reduces  system  capacity. 
Also,  this  method  requires  that  all  transmitters  and  receivers  in  the  system  have  precisely 
synchronized  clocks,  since  just  1  microsecond  of  timing  error  could  result  in  a  300  meter 
position  location  error  [9].  The  use  of  the  TA,  which  is  in  essence  a  timestamp,  is  a 
burden  that  most  would  rather  not  have  to  deal  with.  Furthermore,  the  employment  of  it 
presents  another  problem.  Under  GSM  specifications,  the  TA  is  reported  in  units  of  a  bit 
period,  which  equates  to  a  locus  accuracy  of  554  meters.  Even  this  considerable  amount 
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of  area  will  increase  when  multipath  degradations  and  dilution  of  precision  (the  amount 
by  which  errors  are  degraded  by  geometry)  are  taken  into  account. 

3.  Time  Difference  of  Arrival 

The  TDOA  is  basically  an  improvement  upon  the  TOA  method.  The  idea  behind 
TDOA  is  to  determine  the  relative  position  of  the  mobile  transmitter  by  examining  the 
difference  in  time  at  which  the  signal  arrives  at  multiple  base  station  receivers,  rather 
than  the  absolute  arrival  time.  Therefore,  each  TDOA  measurement  determines  that  the 
transmitter  must  lie  on  a  hyperboloid  with  a  constant  range  difference  between  the  two 
receivers.  The  equation  for  this  range  difference  is  given  by 

R,J  =  -  W  +  O', •  - y?  +  (Z,  - z)2  -  fiTj  - *)!  +  Wj  - yf  +  (Zj  -  z)2  <21) 

where  the  coordinates  (X,  Yj,  Z.)  and  (X,Y.,  Z.)  represent  the  position-fixed  receivers  i 
and  j,  and  x,  y,  and  z  represent  the  unknown  coordinate  of  the  target  transmitter.  [9] 

The  location  of  the  mobile  can  be  estimated  in  two  dimensions  from  the 
intersection  of  two  or  more  hyperboloids  generated  using  TDOA  measurements  from 
three  or  more  base  stations.  This  can  be  extended  to  three  dimensions  if  four  or  more 
TDOA  measurements  are  used. 

As  in  the  TOA  method,  the  TDOA  requires  that  all  fixed  transmitters  and 
receivers  in  the  system  have  precisely  synchronized  clocks.  This,  however,  corresponds 
to  the  timing  standards  already  in  place  at  cellular  base  station  sites.  Unlike  the  TOA 
measurements,  TDOA  measurements  do  not  require  a  timestamp  of  the  transmitted 
signal. 
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III.  POSITION  LOCALIZATION  USING  TDOA  ESTIMATION 


The  process  of  localizing  the  position  of  a  cellular  transmission  occurs  in  two 
stages.  First,  an  accurate  estimate  of  the  TDOA  measured  values  must  be  determined 
from  noisy  signals  that  are  received  from  the  mobile.  This  involves  some  type  of 
denoising  technique,  the  most  popular  being  wavelet  denoising,  which  is  discussed 
thoroughly  in  [10].  Generally  speaking,  this  involves  taking  the  wavelet  transform  of  the 
received  signal,  thresholding  the  wavelet  coefficients,  and  performing  the  inverse  wavelet 
transform  on  the  modified  coefficients.  Second,  a  location  of  position  estimate  is 
determined  from  the  TDOA  calculations.  Both  phases  are  very  complex  and  subject  to 
minor  inconsistencies  that  can  propagate  into  large,  unacceptable  errors.  Thus,  accurate 
and  efficient  algorithms  are  required  for  both  stages  of  processing. 

A.  TDOA  COMPUTATIONS 

In  order  to  be  able  to  evaluate  the  hyperbolic  range  equations,  Eq.  (2.1),  one  must 
first  estimate  the  range  differences  Ry,  or  equivalently,  the  TDOA  tj  -  tj.  The  most 
widely  accepted  method  for  obtaining  these  values  is  the  generalized  cross-correlation 
method.  If  the  transmitted  signal  is  s(t),  then  the  signal  that  arrives  at  one  receiver,  x(t), 
will  be  delayed  and  corrupted  by  noise  and  is  given  by  the  equation 

x(t)  =  as(t-  di)  +  n,(t),  0  <  a  <  1  (3.1) 
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where  dj  is  the  delay  time  from  the  mobile  to  receiver  one,  n;(t)  is  noise,  and  a  is  the  gain. 
Similarly,  the  signal  that  arrives  at  another  receiver,  y(t),  is  given  by  the  equation 

y(t)  =  Ps(t  -  dj)  +  nj(t),  0  <  p  <  1  (3.2) 

where  dj  is  the  delay  time  from  the  mobile  to  receiver  two,  nj(t)  is  the  noise,  and  p  is  the 
gain.  The  cross-correlation  function  between  the  two  received  signals  is  derived  by 
integrating  the  lag  product  of  the  two  signals  for  a  sufficiently  long  time  period,  T.  This 
function  is  given  by 

1  7 

K ,,(0  =  —  J x{t)y(t  -  z)dt  [9],  (3.3) 

This  approach  imposes  the  requirement  that  the  receiving  base  stations  share  a  precise 
time  reference  and  reference  signals,  but  does  not  impose  any  requirement  on  the  signal 
transmitted  by  the  mobile.  We  can  improve  the  SNR  of  the  TDOA  estimation  by 
increasing  the  integration  interval,  T.  The  maximum  likelihood  estimate  of  the  TDOA  is 
obtained  from  the  computed  cross-correlation  function.  This  estimate  is  equal  to  the 
value  of  t  that  maximizes  the  cross-correlation  function,  Eq.  (3.3). 

Another  method  of  estimating  the  cross  correlation  function,  which  yields  the 
same  results,  is  to  compute  the  estimated  cross-spectral  density  function.  An  inverse 
Fourier  Transform  of  the  cross-spectral  density  yields  the  estimated  cross-correlation 
function.  [9] 
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Prior  to  the  cross-correlation  function  computation,  some  frequency  domain 
processing  may  be  used  to  improve  the  ability  to  distinguish  the  desired  signal  from 
arriving  multipath  components.  This  benefit  is  a  result  of  the  conduciveness  of  frequency 
domain  processing  to  filtering  techniques  which  are  effective  in  reducing  the  effects  of 
noise  and  interference  on  the  TDOA  estimate. 

B.  OBTAINING  A  POSITON  LOCATION  ESTIMATE 

After  an  accurate  estimate  of  the  TDOA  measurement  is  obtained,  it  is  then 
possible  to  determine  the  corresponding  solution  to  the  hyperbolic  equation.  This  process 
is  initiated  by  substituting  the  range  difference  estimate  from  Eq.  (2.1)  into  the  hyperbolic 
equation,  Eq.  (3.3),  then  solving  for  the  Cartesian  coordinates  of  the  mobile.  The 
solution  obtained  will  not  be  a  trivial  one,  since  the  range  estimates  will  most  likely  be 
noisy  or  inconsistent.  The  solution  to  the  hyperbolic  equation  is  also  complicated  by  the 
fact  that  the  fixed  location  receivers  are  arranged  in  a  non-uniform  fashion.  This  may  be 
accounted  for  by  linearizing  the  equations  through  Taylor-series  expansion  in  which  only 
the  first  two  terms  are  retained.  However,  this  can  result  in  significant  position  location 
errors  due  to  geometric  dilution  of  precision  (GDOP)  effects.  GDOP  occurs  when  a 
relatively  small  range  of  errors  results  in  a  large  position  location  error  due  to  a  very  large 
difference  in  distance  between  the  mobile  unit  location  on  the  hyperbola  and  the  locations 
of  the  two  receivers  being  used. 
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There  have  been  several  proposed  solutions  to  reduce  the  above  inaccuracies. 
One  being  an  exact  solution  for  the  case  where  the  number  of  equations  equals  the 
number  of  unknown  coordinates.  However,  this  solution  does  not  account  for  redundant 
information.  Other  proposals  include  spherical  interpolation,  a  divide-and-conquer-based 
solution,  and  more  recently,  a  proposal,  by  Y.  T.  Chan,  is  a  closed  form  solution  which  is 
valid  for  both  close  and  distant  sources  [9].  Chan’s  method  has  been  shown  to  be  the 
superior  of  those  mentioned  by  virtue  of  a  lower  computational  complexity  and  higher 
noise  threshold  tolerance.  In  fact,  Chan’s  method  guarantees  optimum  performance 
around  small  TDOA  noise  region  since  it  has  been  shown  to  achieve  the  Cramer-Rao 
lower  bound  in  this  region.  [12] 

The  essence  of  Chan’s  method  is  to  introduce  an  intermediate  variable  to  the 
original  TDOA  equations,  then  transform  the  nonlinear  equations  relating  TDOA 
estimates  and  the  source  position  into  a  set  of  equations  which  are  linear  in  the  unknown 
parameters  and  the  intermediate  variable.  A  least-squares  computation  produces  an 
initial  solution.  By  exploiting  the  known  relation  between  the  intermediate  variable  and 
the  position  coordinates,  a  second  weighted  least  squares  computation  gives  the  final 
solution  for  the  position  coordinates.  [13] 

C.  DETERMINING  THE  ACCURACY  OF  POSTION  LOCATION 

MEASUREMENT 

In  an  effort  to  measure  the  performance  of  TDOA  position  location  systems,  there 
are  several  grading  criterion  at  one’s  disposal,  all  of  which  focus  upon  the  accuracy  of  the 
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position  measurement.  The  most  obvious  and  most  widely  used  measure  of  accuracy  is 
the  root  mean  square  (rms)  accuracy.  This  is  defined  as  the  rms  deviation  of  the 
measured  position  (x,  j>)  about  the  true  position  (x,y),  and  is  given  by  the  following 
equation: 

a  =  ./E[(x-x)2  +(.y-.y)2]  (3-4) 

where  E  is  the  expectation.  The  rms  accuracy  is  a  function  of  the  accuracy  of  the  raw 
locus  measurements  and  the  geometry  of  the  base  stations  relative  to  the  mobile.  There  is 
a  limit  on  the  achievable  accuracy  imposed  by  errors  in  the  locus  measurements.  This 
accuracy  can  be  further  degraded  by  the  physical  geometry  of  the  link  between  base 
stations  to  the  mobile,  referred  to  earlier  as  the  GDOP.  The  GDOP  for  a  2-dimensional 
hyperbolic  position  location  system  is  defined  in  [9]  by: 

I  r 

Jcrx  +  c 

GDOP  =  J-f - y—  (3.5) 

where  g  2  and  a  2  are  the  mean  square  position  location  errors  in  the  x  and  y  directions, 

respectively,  and  aR2  is  the  mean  square  TDOA  ranging  error.  Ideally,  it  would  be 

prudent  to  minimize  GDOP.  However,  to  do  this  would  be  infringent  upon  a  very 

important  principle  in  the  design  of  mobile  communication  systems,  where  the  goal  is  to 
optimize  the  availability  of  service  to  mobile  telephone  users.  Thus,  we  can  minimize  the 
GDOP,  but  we  will  degrade  the  availability  of  service  in  the  process.  Therefore,  the 
focus  of  the  accuracy  determination  must  be  placed  primarily  on  error  in  the  locus 
measurement.  [9] 
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IV.  METHODS  TO  IMPROVE  TDOA  ESTIMATION 


Any  effort  to  improve  the  TDOA  estimation  must  focus  primarily  on  reducing  the 
effective  noise  at  the  receivers.  Earlier  research  presented  seven  different  denoising 
schemes  designed  toward  this  end:  1)  wavelet  denoising  based  on  Donoho’s  method,  2) 
wavelet  denoising  using  the  Wo-So-Ching  threshold,  3)  wavelet  denoising  using 
hyperbolic  shrinkage,  4)  wavelet  denoising  using  median  filtering,  5)  modified 
approximate  maximum-likelihood  delay  estimation,  6)  denoising  based  on  the  fourth 
order  moment,  and  7)  a  time-varying  approximate  maximum-likelihood  technique.  The 
results  of  simulation  of  each  of  these  methods  is  shown  in  Figure  4.1.  [10] 


Figure  4. 1 :  Comparison  of  Denoising  Methods  for  the  GSM  Signal 
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It  was  shown  that  method  1  failed  at  low  SNR’s.  Experimental  results  from  [10] 
showed  that  method  2,  while  an  improvement  over  method  1,  looses  its  advantage  with 
increasing  carrier  frequency.  Additional  improvements  are  made  using  methods  3  and  4, 
but  as  can  be  seen  in  Figure  4.1,  the  best  results  (lowest  mean-square  error)  are  obtained 
using  methods  5  through  7.  We  will  take  a  closer  look  at  these  last  three  methods. 

A.  FOURTH-ORDER  MOMENT  WAVELET  DENOISING 

Inevitably,  any  transmitted  signal  will  acquire  some  type  of  additive  noise  before 
reaching  the  receiver.  However,  it  is  possible  to  estimate  the  frequencies  of  the  corrupted 
complex  sinusoidal  signal.  A  Gaussian  signal  being  completely  characterized  by  its 
second  order  statistics  and  the  odd  order  moments  being  equal  to  zero  for  a  symmetric 
probability  density  function,  the  separation  of  the  signal  and  the  noise  requires  at  least  the 
use  of  the  fourth-order  moments  [11]. 

To  define  the  fourth-order  moment,  we  first  model  the  received  signal  as: 

z(n)  =  x(n)  +  l(n )  (4.1) 

where  x(n)is  a  complex  zero-mean,  non-Gaussian,  fourth  order  stationary  signal  and 
/(«)  is  the  noise,  which  is  a  complex  zero  mean-Gaussian  signal  independent  of  x(n) . 
The  fourth-order  moment  is  then  [1 1]: 

MA:  (j -i,k-  i,l  -  0  ■  E(z,  ,Zj,zk\z')  (4.2) 
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It  was  shown  in  earlier  research  that  the  fourth  order  moment  of  a  detail  function 
which  contains  the  signal  should  be  greater  than  3crnrf  4 ,  where  a nd  4  is  the  noise  power  at 
subband  dj .  [10] 

By  using  this  property,  the  wavelet  coefficients  that  represent  noise  can  be 
eliminated  while  those  having  a  signal  dependency  are  retained.  After  modifying  the 
detail  functions,  the  denoised  signal  is  obtained  by  performing  an  inverse  wavelet 
transform  using  the  modified  coefficients. 

B.  MODIFIED  APPROXIMATE  MAXIMUM  LIKELIHOOD  DELAY 
ESTIMATION 

Critical  to  the  task  of  source  localization  is  the  time  delay  estimation  between 
signals  received  at  two  spatially  separated  sensors  in  the  presence  of  noise.  If  we  let  s(n) 
represent  the  source  signal,  ni(n)  and  n2(n)  represent  the  additive  noises  at  the  respective 
sensors,  and  D  is  the  difference  in  arrival  times  at  the  two  receivers,  then  the  receiver 
outputs,  ri(n)  and  r2(n),  are  estimated  by: 

ri(n)  =  ais(n)  +  m(n),  n  =  0,1,...,T-1  (4.3) 

r2(n)  =  a2s(n  -  D)  +  n2(n),  0  <  ai,  012  <  1  (4.4) 

where  T  is  the  number  of  samples  collected  at  each  channel. 

In  the  approximate  maximum  likelihood  delay  estimation,  after  wavelet 
decomposition  and  prior  to  cross-correlation,  both  the  channel  outputs  are  optimally 
weighted  at  different  frequency  bands  with  the  use  of  an  orthogonal  wavelet  transform. 
This  weighting  is  done  to  attenuate  the  noise  spectral  components.  The  scaled  subband 
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components  are  then  combined  using  inverse  wavelet  transform  to  construct  the  Modified 
AML  (MAML)  prefiltered  signal.  The  orthogonal  wavelet  transform  is  attractive  because 
in  addition  to  being  computationally  efficient,  it  does  not  suffer  from  the  performance 
degradation  due  to  errors  that  may  occur  in  the  estimation  of  D  used  in  Eq.  (4.4),  since  in 
practice,  a  highly  accurate  estimate  of  D  may  be  difficult  to  achieve.  [12] 

The  final  MAML  delay  estimate  is  calculated  from  the  location  of  the  peak  of  the 
cross-correlation  function  of  the  two  denoised  signals. 

C.  TIME-VARYING  MODIFIED  APPROXIMATE  MAXIMUM 

LIKELIHOOD  DELAY  ESTIMATION 

The  time-varying  MAML  method,  as  is  obvious  from  the  name,  works  in  much 
the  same  way  as  the  MAML,  except  that  it  additionally  accounts  for  the  time-varying 
characteristics  of  the  signal.  Due  to  the  time-varying  nature  of  some  signals,  the 
subbands  may  not  contain  actual  signal  components  throughout  the  entire  segment. 
Therefore,  the  MAML  method  may  loose  some  of  its  effectiveness  when  applied  to 
segments  in  the  presence  of  transitory  signals.  The  time-varying  MAML  method,  in 
essence,  allows  for  the  MAML  method  to  be  applied  twice  to  each  segment,  thus 
approximating  a  time-varying  gain. 
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V.  SIGNAL  SIMULATION  AND  TEST  RESULTS 


A.  SIGNAL  GENERATION 

All  signals  used  in  this  project  for  test  purposes  were  generated  using  the  Hewlett- 
Packard  Advanced  Design  System  (HP-ADS).  This  provided  an  excellent  method  of 
obtaining  signals  that  are  essentially  the  same  as  would  be  encountered  in  an  actual 
cellular  communication  receiver. 

1.  Advanced  Design  System 

The  HP-ADS  program  is  an  invaluable  tool  for  the  engineers  representing  many 
different  design  aspects  such  as  communications,  digital  signal  processing,  electronic 
circuits,  mechanical  circuits,  and  many  others.  Of  course,  for  this  project,  we  were 
particularly  interested  in  the  communications  package.  This  package  allows  the  user  to 
custom  design  any  type  of  communications  system,  run  simulations  of  the  design,  and 
extract  data  collected  from  the  simulation.  Figure  5.1  shows  the  system  used  for  the 
signal  simulation.  A  more  detailed  view  of  each  of  the  major  components  of  this  system 
can  be  found  in  Appendix  [B].  Appendix  [B]  also  provides  detailed  instructions  for  using 
the  HP-ADS  program  to  design  and  simulate  a  communications  system. 

2.  GSM  Signal  Generation 

One  parameter  that  must  be  specified  for  the  GSM  signal  generation  is  the  sample 
time.  The  specifications  of  the  GSM  signal  were  presented  in  Chapter  II.  The  HP-ADS 
program  specifies  a  symbol  period  of  3.7037  microseconds.  The  filter  bandwidth  is  1.2 
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MHz.  This  allows  us  to  sample  at  a  minimum  frequency  of  2.4  MHz  without  violating 
the  Nyquist  criterion.  For  test  purposes,  we  use  three  different  sampling  frequencies.  As 
shown  in  Table  5.1,  these  are:  10.8  MHz,  which  corresponds  to  a  sample  time  of  92.592 
microseconds;  5.4  MHz,  which  corresponds  to  a  sample  time  of  185.185  microseconds; 
and  2.7  MHz,  which  corresponds  to  a  sample  time  of  370.370  microseconds. 


Figure  5.1:  HP-ADS  GSM  Communications  System 
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Samples/Symbol 

Sample  Time 
(microseconds) 

Sampling  Frequency 
(Megahertz) 

Symbols/600 

Samples 

40 

92.592 

10.8 

15 

20 

185.185 

5.4 

30 

10 

370.37 

2.7 

60 

Table  5.1:  GSM  Signal  Parameter  Combinations 

The  I  and  Q  channel  outputs  of  each  of  these  GSM  signals  are  shown  in  Figure 
5.2.  One  fact  that  is  immediately  evident  from  Figure  5.2  is  that  as  the  sampling  interval 
increases,  so  does  the  variation  in  the  GSM  signals.  The  reason  for  this  is  obvious  if  we 
refer  to  Table  5.1.  We  can  see  that  as  the  sample  interval  increases,  we  can  expect  to 
contain  more  symbols  within  a  given  number  of  samples.  For  test  purposes,  we  use  a 
constant  number  of  samples  (600)  for  all  simulations.  We  also  realize  a  better  correlation 
characteristic  with  increasing  sample  time  as  shown  in  Figure  5.3.  Better  correlation  is 
demonstrated  by  a  sharper  main  lobe  with  lower  sidelobes.  Figure  5.3  (d)  shows  the 
correlation  of  the  Matlab  generated  GSM  signal  presented  by  Aktas  [10].  This  figure 
represents  the  best  correlation.  Plots  of  the  power  spectral  densities  of  these  signals  can 
be  found  in  Appendix  [C]. 
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(b)  I  channel;  Ts  =  92.592  ns 


(a)  Q  charnel;  Ts  =  92.592  ns 


(b)  Corr^ationof  sigTaisxandy;  ADS185.185rsec 
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Figure  5.3:  Auto-Correlation  Results  of  Test  Signals 

B.  SIMULATIONS 

In  Chapter  IV,  we  discussed  three  methods  to  improve  the  TDOA  estimate.  In 
this  section,  we  will  use  these  methods  to  test  our  HP- ADS  generated  signals  and  see  if 
we  are  able  to  extract  useful  TDOA  data.  Bear  in  mind  that  our  criterion  for 


(c)  Correlation  of  signals  x  and  y;  ADS  37Q.37nsec 


(d)  Correlation  of  signals  x  and  y;  MAUAB 


(a)  Correlation  of  signals  x  and  y;  ADS  92.59  nsec 


improvement  in  TDOA  estimation  is  a  lower  mean-squared  error.  Each  method  was  tried 
using  four  different  realizations  for  each  of  the  three  sampling  times,  for  a  total  of  twelve 
different  realizations  per  method  used.  However,  we  will  only  present  one  set  from  each 
sample  time  here,  as  each  set  follows  a  similar  trend.  The  remainder  of  the  plots  can  be 
found  in  Appendix  [D].  To  simplify  explanation,  we  shall  label  the  data  sets  as  follows: 

a)  HP- ADS  generated  signal  with  a  sample  interval  of  92.59  nsec. 

b)  HP- ADS  generated  signal  with  a  sample  interval  of  1 85. 1 85  nsec. 

c)  HP-ADS  generated  signal  with  a  sample  interval  of  370.370  nsec. 

d)  Matlab  generated  signal. 

We  will  use  the  Matlab  generated  signal  as  our  benchmark  for  comparison.  For  each 
method  we  will  plot  the  mean-squared  error  against  SNR’s  in  the  range  of  -6  dB  to  20 
dB.  All  Matlab  codes  used  in  these  simulations  are  provided  in  Appendix  [E].  For 
clarity  in  the  following  discussion,  we  define  the  terms  average  error  and  percentage 
improvement.  Average  error  is  calculated  by  adding  all  non-zero  error  values  for  a 
particular  realization,  and  dividing  that  sum  by  the  total  number  of  the  non-zero  elements 
in  that  realization.  Percentage  improvement,  as  shown  in  the  equation  below,  is 
calculated  by  subtracting  the  improved  value  from  the  original  value  and  dividing  the 
difference  by  the  original  value  then  multiplying  by  100. 

original  value  -  improved  value 

— — - — - x  100%  =  percent  improvement 

original  value 
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mean  square  error  m ean  square  error 


1 


Fourth  Order  Moment 


The  results  using  this  method  are  shown  in  Figure  5.4.  There  are  two  obvious 
trends  that  are  noticed  in  this  figure.  We  can  easily  see  that  the  mean-squared  error 

(a)  4th  Order  Moment  HP-ADS;  Ts  =  92.59  nsec  (b)  4th  Order  Moment  HP-ADS;  Ts  =  185.185  nsec 
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(d)  4th  Order  Moment;  Matlab 


xcorr  without  denoising 
4th  order  moment 


0  5  10  15 

SNR  (dB) 


Figure  5.4:  Results  of  Varying  Sampling  Interval  for  Fourth-Order  Moment  Technique 


improves  (decreases)  as  the  sampling  interval  increases.  This  observation  is  presented 
numerically  in  Table  5.2.  From  this  table,  we  can  clearly  see  a  decrease  in  error  values  as 
we  move  from  left  to  right  across  each  row. 

In  the  Matlab  plot  (Figure  5.4  (d)),  or  our  benchmark  plot,  notice  that  the  mean- 
squared  error  is  zero  at  6  dB  and  higher  SNR  for  the  Fourth-Order  Moment  denoised 
curve.  Also  notice  that,  as  expected,  the  denoised  curve  is  lower  than  the  non-denoised 
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2.66 

1.34 

0.70 

0.31 

0.08 
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9.03 

3.82 

2.11 
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0.88 

0.35 

0.11 

0.02 

6 

2.89 
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0.67 

0.49 

0.10 

0 

0 

9 

1.11 

0.86 

0.65 

0.27 

0.19 

0.01 

0 

0 

12 

0.63 

0.24 

0.240 

0.05 

0.02 

0 

0 

0 

15 

0.17 

0.05 

0.060 

0 

0 

0 

0 

0 

18 

0.03 

0.01 

0 

0 

0 

0 

0 

0 

20 

0 

0 

0 

0 

0 

0 

0 
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Table  5.2:  Mean-Squared  Error  Values  at  Varying  SNR  for  the  Fourth  Order  Moment 
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curve  at  all  points.  This  shows  that  the  denoising  is  effective  in  producing  an  improved 
TDOA  estimate. 

Using  the  HP-ADS  data  with  a  sampling  interval  of  92.59  nsec,  case(a),  the  mean- 
square  error  at  -3  dB  SNR  is  26.3  when  fourth-order  denoising  is  employed.  This  is  a 
3 1 .4%  improvement  from  the  non-denoised  curve.  The  average  percentage  improvement 
over  all  SNR’s  between  -3  and  20  is  45.6%. 

Using  HP-ADS  data  with  a  sampling  interval  of  185.185  nsec,  we  notice  an 
average  52.2%  increase  with  fourth-order  moment  denoising.  The  mean-squared  error  is 
5.68  at  -3  dB,  which  is  a  78.4%  improvement  over  the  -3  dB  value  of  case  (a). 

Using  HP-ADS  data  with  a  sampling  interval  of  370.37  nsec,  case  (c),  we  obtain 
even  better  results.  The  mean-squared  error  at  -3  dB  is  now  1.96,  which  is  a  65.5% 
improvement  over  case  (b).  The  average  improvement  realized  by  using  fourth-order 
moment  denoising  instead  of  the  non-denoised  method  in  this  case  is  69.0%.  We  now 
see  a  definite  trend  that  as  our  sampling  interval  increases,  denoising  increasingly 
improves  the  mean  squared  error  values. 

In  the  benchmark  case  (d),  there  is  a  maximum  mean-squared  error  of  0.25  at  -3 
dB  using  fourth-order  denoising.  The  average  percentage  improvement  of  fourth-order 
moment  denoising  over  non-denoising  is  72.5%.  We  also  notice  that  the  mean-squared 
error  is  zero  for  all  SNR  above  6  dB. 

The  improvements  we  have  discussed  in  this  section  can  be  related  to  the 
correlation  function  plots  presented  in  Section  5. A.  There,  we  saw  that  the  correlation 
function  plot  improved  with  increasing  sample  time.  Here,  we  state  that  a  higher  degree 
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mean  square  error  mean  square  error 


of  correlation  of  the  signal  components  reduces  the  probability  of  error  in  the  TDOA 
calculations.  Thus,  we  achieve  lower  mean  square  error  values  with  increasing  degree  of 
correlation  of  signal  components. 


2.  Modified  Approximate  Maximum  Likelihood  (MAML) 

The  results  using  this  method  are  shown  in  Figure  5.5.  We  notice  a  similar  trend 


as  in  the  fourth  order  moment  method.  We  end  up  with  lower  mean  squared  error  values 


(a)  MAML  HP-ADS;  Ts  =  92.59  nsec  (b)  MAML  HP-ADS;  Ts  =  185.185  nsec 


SNR  (dB)  SNR  (dB) 


SNR  (dB)  SNR  (dB) 

Figure  5.5:  Results  of  Varying  Sample  Times  for  MAML  Technique 
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as  the  sampling  interval  increases.  We  also  notice  that  these  mean  squared  error  values 
are  lower  than  those  obtained  in  the  fourth  order  moment  method,  as  evidenced  in  Table 
5.3. 


In  caSe  (a),  there  is  an  average  improvement  in  mean  squared  error  of  89.6%  by 
employing  MAML  denoising.  Recall,  from  Section  5.B.1,  the  fourth  order  moment 
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Table  5.2:  Mean-Squared  Error  Values  at  Varying  SNR  for  MAML  Method 
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denoising  method  produced  only  an  average  45.6%  improvement  for  case  (a).  Thus,  we 
have  proved  that  the  MAML  denoising  method  is,  in  fact,  superior  to  the  fourth  order 
moment. 

In  cases  (b)  and  (c),  there  is  an  average  improvement  of  89.5%  and  88.8% 
respectively.  These  are  in  accord  with  the  averages  obtained  in  case  (a).  Thus,  we  have 
shown  that  through  the  use  of  the  HP-ADS  program,  we  can  generate  a  simulated  GSM 
signal  and  extract  the  necessary  data  to  obtain  a  TDOA  estimation. 

3.  Time-Varying  Approximate  Maximum-Likelihood  (TV-MAML) 

In  Figure  4.1,  we  saw  that  the  time-varying  MAML  (TV-MAML)  method 
produced  the  best  results  at  higher  SNR.  The  data  extracted  from  the  HP-ADS  program 
did  not  produce  similar  results  as  obtained  in  the  fourth  order  moment  technique  and  the 
MAML  method.  An  unsuccessful  run  of  the  TV-MAML  program  is  shown  in  Figure  5.6. 
We  can  see  in  this  figure  that  the  MAML  plot  outperforms  the  TV-MAML  plot. 
Furthermore,  the  TV-MAML  plot  has  a  non-zero  mean  squared  error  even  at  very  high 
SNR.  This  can  be  explained  as  follows.  Recall  from  Section  4.C.,  the  TV-MAML 
requires  the  MAML  technique  be  applied  twice  to  each  segment  of  the  signal  in  order  to 
reduce  the  probability  of  capturing  non-signal  components  along  with  the  actual  signal. 
Recall  also,  all  of  our  simulations  used  data  sets  that  were  600  samples  long.  Therefore, 
the  degree  to  which  a  particular  data  set  is  transitory  depends  upon  its  sample  interval. 
Lower  sampling  intervals,  as  described  in  Figure  5.2,  produce  data  sets  that  contain  less 
information.  Thus,  the  potential  effectiveness  of  the  TV-MAML  is  not  fully  realized  if 
the  signal  does  not  contain  a  relatively  large  amount  of  information. 
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TV -AML  Using  HP-ADS  Data  Ts  =  370.37  nsec 


-5  0  5  10  15  20 


SNR 

Figure  5.6:  Results  of  Unsuccessful  Run  of  TV-MAML  Program  Using  Sample  Interval 
of  370.37  nsec 

In  order  to  evaluate  the  effectiveness  of  the  TV-MAML  method,  we  must  use  a 
higher  sampling  interval  when  collecting  data  from  the  signal.  Thus,  we  found  740 
nanoseconds  to  be  a  suitable  sample  time.  The  signal  is  shown  in  Figure  5.7  along  with 
the  results  from  the  TV-MAML  program.  Here,  we  see  a  definite  improvement  relative 
to  Figure  5.6.  The  TV-MAML  outperforms  the  MAML  program  at  0  dB  (and  lower 
SNR’s)  by  an  average  of  18.1%.  This  demonstrates  that  the  TV-MAML  works  better  in  a 
higher  noise  environment  than  the  MAML. 
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VI.  SUMMARY  AND  CONCLUSIONS 


A.  SUMMARY 

The  objectives  of  this  thesis  are  to  simulate  GSM  signals  and  to  evaluate  the  mean 
squared-error  for  TDOA  estimation.  The  Hewlett-Packard  Advanced  Design  System  is  a 
very  powerful  communications  development  tool  and  proved  invaluable  toward  our 
objective  of  GSM  signal  generation.  Appendix  [F]  contains  a  concise  user’s  guide  to  the 
HP-ADS  program  for  performing  the  tasks  conducted  in  this  thesis.  We  were  able  to 
manipulate  this  system  to  provide  signals  with  desired  characteristics  which  would  later 
be  used  in  the  Matlab  environment  to  determine  the  effectiveness  of  three  different 
denoising  methods.  The  denoising  methods  were  presented  in  earlier  research.  These  are 
the  fourth-order  moment,  modified  approximate  maximum-likelihood,  and  time-varying 
approximate  maximum-likelihood  methods.  We  used  three  signal  sets  generated  by  the 
HP-ADS  program  to  test  the  performance  of  each  of  these  methods.  The  results  of  these 
tests  agreed  with  the  earlier  research.  It  showed  that  the  MAML  method  is  the  best  choice 
except  when  the  signal  contains  a  relatively  large  amount  of  information,  in  which  case 
the  TV-MAML  tends  to  be  better. 

In  the  process  of  testing  the  denoising  techniques,  we  have  demonstrated  our 
ability  to  extract  localization  information  from  our  generated  GSM  signals.  This  is 
manifested  in  the  fact  that  we  obtain  relatively  low  mean  squared  error  in  our  TDOA 
calculations.  We  also  found  that  the  error  could  be  reduced  further  by  increasing  the 
interval  over  which  samples  are  taken  from  the  signal. 
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Although  the  Matlab  programs  provided  the  desired  results,  there  were  some 
inconsistencies  to  be  noted.  In  the  Matlab  code  for  the  TV-MAML,  wavelet 
decomposition  and  re-composition  is  performed  on  the  noise-infected  input  signal.  At 
some  point  between  the  decomposition  and  re-composition,  a  dc  offset  is  introduced  as 
well  as  some  slow  frequency  oscillation  toward  the  end  of  the  sequence.  These  issues 
will  be  addressed  in  subsequent  work  on  this  project. 

B.  FUTURE  WORK 

In  addition  to  providing  a  solution  to  the  oscillations  and  DC  offset  described  in 
Section  6.A.,  it  may  be  possible  to  modify  the  existing  Matlab  code  to  obtain  even  better 
mean  squared  error  values.  Follow  on  work  should  also  re-examine  the  limitations  of  the 
time- varying  approximate  maximimum-likelihood  method  presented  in  Section  5.3.  An 
extension  to  cases  with  different  SNR’s  is  also  needed  to  determine  the  performance  in  a 
fading  environment.  Also,  the  situation  where  there  is  a  weak  signal  in  one  channel  and  a 
stronger  signal  in  the  other  channel  should  be  examined.  This  would  provide  insight  into 
the  question  of  what  ranges  of  SNR  are  reliable.  Finally,  an  investigation  into  the  actual 
position  location  estimation  using  the  solution  of  the  hyperbolic  equation  presented  in 
Chapter  HI  should  be  conducted. 
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APPENDIX  A.  USA  TODAY  ARTICLE  ON  CELLULAR  911  LEGISLATION 


10/26/99-  Updated  05:39  PM  ET 

Clinton  signs  nationwide  911  law 

WASHINGTON  (AP)  -  President  Clinton  signed  legislation  Tuesday 
making  91 1  the  official  emergency  number  nationwide  -  for  both 
regular  and  cellular  phones. 

The  measure  also  cadis  for  development  of  technology  that  can  track 
mobile  callers. 

People  with  wireless  phones  now  will  be  able  to  speed  responses  to 
highway  accidents,  crimes  and  natural  disasters,"  Clinton  said. 
"Getting  rapid  care  tb  someone  who  is  suffering  from  a  heart  attack  or 
is  involved  in  a  car  crash  can  mean  the  difference  between  life  and 
death." 

While  91 1  is  widely  used  as  the  emergency  number  for  traditional 
phones,  there  are  20  different  codes  for  wireless  callers  across  the 
country.  The  changes  are  aimed  at  cutting  response  times  for  the 
crews  who  answer  98,000  emergency  calls  daily  from  cellular  phone 
callers. 

"In  my  home  state,"  said  Sen.  Conrad  Bums,  R-Mont,  "three  quarters 
of  the  deaths  in  rural areas  are  because  the  first  responders  couldn’t 
get  there  in  time." 

Health  care  professionals  joined  Bums  at  a  Capitol  Hill  news 
conference  to  applaud  the  new  law. 

"We  have  great  emergency  room  personnel.  We  can  do  a  lot  for 
,  accident  victims  if  we  can  find  them  and  gefthem  there,"  said 
Barbara  Foley  of  the  Emergency  Nurses  Association.  "That’s  what 
this  legislation  helps  us  do." 

Another  provision  of  the  act  directed  the  Federal  Communications 
Commission  to  help  states  develop  emergency  systems,  including 
technology  that  can  automatically  locate  cellular  callers  who  have 
dialed  91 1  or  been  involved  in  an  accident 

The  FCC  in  September  moved  forward  with  plans  to  require  that 
cellular  91 1  calls  automatically  provide  a  caller's  location.  Regulators 
want  manufacturers  to  begin  providing  locator  technology  within  two 
years. 
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Privacy  advocates  have  raised  concerns  about  potential  abuse  of  the 
technology,  which  would  take  advantage  of  the  Global  Positioning 
System  developed  by  the  military. 

The  law  signed  Tuesday  called  on  regulators  to  establish  "appropriate 
privacy  protection  for  call  location  information,"  including  systems 
that  provide  automatic  notification  when  a  vehicle  is  involved  in  an 
accident. 

It  said  that  calls  could  only  be  tracked  in  nonemergency  situations  if 
the  subscriber  had  provided  written  approval.  "The  customer  must 
grant  such  authority  expressly  in  advance  of  such  use,  disclosure  or 
access,"  according  to  Senate  documents  detailing  provisions  of  the 
legislation. 

An  estimated  700  small  and  rural  counties  have  no  coordinated 
emergency  service  to  call  -  even  with  traditional  phones.  The  bill 
would  encourage  private  91 1  providers  to  move  into  those  areas  by 
granting  the  same  liability  protections  to  wireless  operations  that  now 
are  offered  to  wireline  emergency  service  systems. 

Separately,  the  FCC  took  action  earlier  this  year  to  increase  the 
number  of  cellular  calls  to  91 1  that  are  successfully  completed.  The 
commission  required  that  new  analog  cellular  phones  -  not  existing 
phones  -  be  made  with  software  that  routes  91 1  calls  to  another 
carrier  when  a  customer's  own  service  cannot  complete  the  call. 

Calls  sometimes  aren't  completed  because  a  caller  is  in  an  area  where 
his  or  her  carrier  does  not  have  an  antenna,  because  networks  are 
overloaded  or  because  buildings  or  geography  block  signals. 

Digital  phones,  of  which  18.8  million  now  are  in  use,  were  not 
covered  by  the  new  FCC  rules  adopted  in  May  because  such  phones 
are  more  complex  than  their  analog  counterparts  and  there  is  no  easy 
fix  for  the  problem. 


•  Go  to  Washington  news 

•  Go  to  News  frontpage 
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APPENDIX  B.  HP-ADS  GSM  COMMUNICATION  SYSTEM  COMPONENTS 


MODULATOR 


TRANSMITTER 


NJones 

XI 

TStefFTStep 
Frequency  1*330  MHz 

Phase  1=0.0 

Addilki<alTcn«s**r 

Randont®nase=Ko 

PnaseNciseDa  la- ^ 

«OaP50.OOlim 

RTemp=-274.C- 

?N-T>pe=Randfifn  PN 


czr~z 

!K..' 

jj;sj 

»  p<  pi ;  j  i 

!j~i 

j  ji^j 

EPFJ&aeiwoftliTimed 
61  “ 

Rin=50.0  Ohm 

R0«1=$3.00fim 

RTeinp=-27*,0 

FCcrt1er*900  WKz 

Pw55andwK«h=1^  MHz 

PassAt&isO.C 

SSopBafidvrtfflft*4MHz 

SlopAttejisSO.O 

N*3 

.VmtfewTypfi^Rcaangular 
irtpTvn«=C-0  jec 


GainR? 

Gi 

R^SO.OOhm 
ROyi=5C.O  Ohni 

Gain=b+/t.O 

Ncnseftg^U 

GC7ype=70i+tJ3i1 

7Clwt=&nic*{-l3} 

P$at=1.0W 
GCSal'1.0 
GComp='*C.O  0.D  0.0  * 


Rm=50-C  Ohm 

ROut=£G.O  Ohm 

RTemp=-27<.0 

G?kv=22+j'0.0 

Nct$eFigun*3.5 

GC7ype=TOM3c1 

70Hu1*dDmto<-3) 

dBclou&i&miwt-iO) 

PSat*1.0W 

GCSa^l.O 

GCoflip=’0M00.0* 


MixerR? 

Ml 

Sin=50.0Ghm 

RO^=50.0  Ohm 

RTmp=*2?4.3 

7ype*RFmint>sLC 

RfRe^KO.O 

lmRej-203,0 

lcRij^200.0 

N-oiseFigurB*? 

lCamp=‘-6,r 


SPFJutterwc&fiTimed 

32 

RIr.^50.0  Ohm 

R04=KXIOhm 

RTernpJS-274,!> 

FCcafertCMHz 

PassB3n6«Hh=1.2MHz 

PassA3en*3;0 

sia?B3ftewdm=<  MHz 

StopAaen^C 

N=3 

WindcwType^Reciansular 

ImpTim^O.OiK; 


GainR  F 
G3 

Rii^so.o  onm 

ROi#:$a.OOhm 

RTemp=‘274.0 

Gni?F23+j*0,0 

NoJsengum=3i 

GCTypesTOWBcl 

70kw!=dbmtow{5> 

d8cloi!t=<JKntow(*5] 

PSat=1.CW 

GC$35=1.0 

GCom?=U,00,0  0.0’ 


GairflF 

G< 

Rft=50.0  Ohm 

ROuv-SO.O  Ohm 

RTfmF-274.0 

Gaifi=2$*j‘0.3 

NeissFiguns^ 

GCType'noiie 

TO&1=d&nno#(5) 

d3cloitafcmtc-*(*5j 

P$35*1.0W 

GC$ai=1.0 

GCotp='0.0  0,0  Q.C  * 
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RECEIVER 


CONN; 

Num=1 

layer-cond 


O- 

Port 

C0NN1 

NuitfI 

layer=cond 


_ ,'N  Tones 

\;  InT 

\  !  TStep=TStep 

_ I  Frequency-)  =830  MHz 

Power1=5  mW 
w  Phasei=0.0 

Addition  alTone$= 
Random  Phase= No 
PhaseNoiseData= 

'f  R0ut=50.0  Ohm 


RTemp=-274.0 
PN  Tvoe-Random  PN 


MixerRF 

M2 

Rln=50.0  Ohm 
R0ut=50.0  Ohm 
RTemp=-274.0 
Type=RF  plus  LO 
RfRej=-200.0 
lmRej=-200.0 
LoRej=-200.0 
Noi$eFjgure=8 
LComp—  -7 , 5 


BPF  ButterworthTimed 
B1  “ 

Rln=50.0  Ohm 
R0ut=50.0  Ohm 
RTemp=-274.0 
FCenter=S00  MHz 
PassBandwidth='!.2  MHz 
PassAtten=3.0 
StcpBandwidth=80  MHz 
StopAtten=50.0 
N=3 

WindowType=Rectangular 

lmpTime= 


i - ; 

r-— 1 

!;>v  ? 

'  > 

— ►* — 

I. 

y 

V 

-o 

~  ^  Pen 

CONN2 

GainRF 

GainRF 

Num=2 

Cl 

G2 

laver-coi 

Rln=50.00hm 

Rin=50.0  Ohm 

ROut=50.0  Ohm 

ROut=50.0  Ohm 

RTemp=-274.0 

RTemp=*274.0 

Gain=3+j’0.0 

Gain=25+j*0.0 

Noise  Figure=3.5 

NoiseFigure=6 

GCType=none 

GCType=TOI*dBd 

T0lout=3.0  W 

TOIout=dbmtow(37) 

dBc1cut=1.0  W 

dBc1out=dbmtow(27) 

PSat=1,0W 

P$at=l.0W 

GCSat=1.0 

GCSat=i.O 

GCorr.p="0.0  0.0  0.0 " 

GComp="0.0  0.0  0.0 " 

DEMODULATOR 


BPF^RaisedCosineTimed 
B1  ~ 

Rln=50.0  0hm 
ROut=50.0  Ohm 
RTemp=-274.0 
FCenter=69.7  MHz 
Bandwidth=1  MHz 
Exce$sBw=0.0 
Type=lmpulse  model 
SquareRoot=No 
Delay=10  usee 
WindowType=Rectanguiar 


SplitterRF 

SI 

Rln=50.0  Ohm 
ROut=50.0  Ohm 
RTemp=-274.0 


GMSKJRecovery 

G1 

Rin=50.0  Ohm 
R0ut=50.0  Ohm 
RTemp=-274.0 
RefFreq=70  MHz 
SymbolTime=Symbornme 
RecoveryBw=0.01  MHz 
- 


TT 


"T5  cT< 

DEM 

3MSK 


GMSK_Demod 

G2 

Rln=50.0  Ohm 
ROut=50.0  Ohm 


RTemp=-274.0 


<3 

Port 

C0NN2 

Num=2 

layer=cond 
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APPENDIX  D.  MATLAB  SIMULATIONS  FOR  EACH  DATA  SET 


FOURTH  ORDER  MOMENT  &  MAML  FOR  SAMPLE  TIME  =  92.592  nsec 


Results  for  DATA  92  setl 


Results  for  DATA  92  set3 


1  .  -1 . -  ■  1 

—  xcorr  without  denoising 
-  4th  order  moment 

—  Modified  AML 

0  5  10  15  2( 


Results  for  DATA  92  set2 


Results  for  DATA  92  set4 


SNR  in  (dB) 


SNR  in  (dB) 


FOURTH  ORDER  MOMENT  &  MAML  FOR  SAMPLE  TIME  =  185.185  nsec 
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APPENDIX  E.  MATLAB  CODES 


FOURTH-ORDER  MOMENT  DENOISING  TECHNIQUE 


************ 

<^********************************************************** 

************ 

%  Denoise_sta:  Wavelet  Denosing  Based  on  The  Fourth  Order 
Moment 
% 

%  SYNTAX:  y=Denoise_sta (xn, yn) 

% 

%  INPUT:  xn  =  Received  signal  from  first  receiver 
%  yn  =  Received  signal  from  second  receiver 

q, 

'o 

%  OUTPUT:  y  =  Denoised  signal  Xn  based  on  Yn  statistics 

o, 

o 

%  SUB_FUNC:  None 
%  Written  by  Spiros  Mantis 

<^********************************************************** 

************ 

•^c'k-^r-krk’k'k'k'k'k'k'k 

function  y=Denoise_sta (xn, yn) ; 

xyn=xcorr (xn, yn, ’biased’ ) ; 

[sigmas  b]=max(xyn); 
rx=xcorr (xn, ’ biased ’ ) ; 
maxx=rx (length (xn) ) ; 
ry=xcorr (yn, ’biased’ )  ; 
maxy=ry (length (yn) ) ; 
s igmanl=maxx- sigmas; 
sigman2=maxy-sigmas; 

lamdax=3 . l*sigmanlA2; 
lamday=3 . l*sigman2A2; 

nx=floor (log2 (length (xn) ) ) ; 
ny=floor (log2 (length (yn) ) ) ; 
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[cx  lx] =wavedec (xn, nx, ' db4 ' )  ; 
[cy  ly] =wavedec (yn, ny, ' db4 ' ) ; 

dxc=  [  ] ; 
for  i=l:nx 

d=detcoef (cx, lx, i) ; 
dl=length (d) ; 

A= ( 1/dl) *sum (d. A4 )  ; 

if  AClamdax 

dc=zeros (1, dl) ; 
else 

dc=d; 

end 

dxc=[dc  dxc] ; 

end 

a=appcoef (cx, lx, ' db4 ' , nx) ; 
al=length (a) ; 

A= (1/al) *sum(a. A4) ; 

if  A<lamdax 

ac=zeros (1, al) ; 
else 

ac=a; 

end 

dxc=[ac  dxc]  ; 

xd=waverec (dxc, lx, ' db4 ' ) ; 

y=xd; 
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2-********************************************************** 

o 

o 

'k'k'k'k'k'k'k'k'k'k'k'k 

%  Modif ied_tez6 :  This  is  a  test  program  for  wavelet 
denosing 

%  based  on  the  fourth  order  moment  technique. 

%  GSM  signal  is  used. 

% 

%  SYNTAX:  Modif ied_tez6 
% 

%  INPUT:  None 
% 

%  OUTPUT:  Mean  square  error  versus  SNR 
% 

%  SUB_FUNC:  Denoise_sta .m 
%  Written  by  Spiros  Mantis 

%**it**i?*i'****ji:********ic***icjc****icjr*******'k'*:'k'k'k-k'*rj<:jc'k'k*ic'k'k'k'k'k 

O 

-k-k'jc-jc'kJcJt'k'k'kJc'k 

o 

•k'k'k'k'kic'k^r'k'k'k'k 

clear  all 

%gsm_set;  %Conf iguration  variables  created  in  memory. 

%  these  are: 

%  Tb (=  3 . 692e-6) 

%  BT(=  0.3) 

%  OSR (=  4) 

%  SEED (=  931316785) 

%  INIT_L(=  260) 

%data=data_gen (INIT_L) ;  %  this  creates  a  binary  data 
% [tx_burst, I, Q]=gsm_mod(Tb, OSR, BT, data, TRAINING) ; 

%s=I+j  *Q; 
data_370_set4 
s=transpose  (s2)  ; 
sl=length (s) ; 

pow= (1/sl) *sum(abs (s) . A2) ; 


K=100  %  number  of  realizations 

rand ( ' seed' ,40) ; 
f=150*rand (1, K) ; 
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delay=floor (f ) ; 

delay (l:K/2)=-delay (l:K/2) ;  %  delay  is  between  -150  to  +150 


n= [20  18  17  16  15  14  13  12  9  6  3  0  -3  -6]  ; 

SNR=10.A (n./lO)  ; 

for  k=l : length (SNR) 

oran=SNR (k) 
for  i=l:K 

x= [zeros (1, 200)  s  zeros  (1, 224 )] ; 

y= [zeros (1, 200+delay(i) )  s  zeros (1, 224-delay (i) )] ; 
randn ( ' state ' ,2* (i+j )  )  ; 

noil_real=sqrt (pow/ (2*oran) ) *randn (1,  1024) ; 
randn ( ' state ' , 3* (i+j )  )  ; 

noil_imag=sqrt (pow / (2*oran) ) *randn (1, 1024) ; 
randn (' state 4* (i+j )) ; 

noi2_real=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 
randn ( ' state ' , 5* (i+j ) ) ; 

noi2_imag=sqrt (pow/ (2*oran) ) *randn(l,  1024) ; 

noil=noil_real+j  *noil_imag; 
noi2=noi2_real+j  *noi2_imag; 

xn=x+noil;  %  x  +  noise 
yn=y+noi2;  %  y  +  noise 


%  TDOA  calculation  with  xcorr (  without  denoising) 

xy=xcorr (xn, yn) ;  %  correlation  of  x  and  y  with  xcorr 
[al  bl]=max(real(xy)); 

erl (i) =delay (i)-(bl-1024) ; 


%  Fourth  order  moment  denoise 


60 


X_real=Denoise_sta  (real (xn) , real (yn) ) ; 
X_imag=Denoise_sta (imag(xn)  ,  imag (yn) ) ; 
Y_real=Denoise_sta (real (yn)  ,  real (xn) ) ; 
Y_imag=Denoise_sta (imag (yn)  , imag(xn) ) ; 

X=X_real+ j  *X_imag; 

Y=Y_real+j  *Y_imag; 

XY=xcorr (X, Y) ;  %  Correlation  of  X  and  Y  (denoised) 

[a2  b2] =max (real (XY) ) ; 

erlO (i) =delay (i) - (b2-1024) ; 


end 

errorl0a370set4 (k) = (1 /length (er 10) ) *sum (erlO . A2) ; 
errorlsta370set4 (k) = (1 /length (erl) ) *sum(erl. A2)  ; 

H10a370set4 (k,  : ) =erl0; 

Hsta370set4 (k,  : ) =erl; 


end 


figure (6) 

k= [20  18  17  16  15  14  13  12  9  6  3  0  -3  -6] ; 

plot (k, errorlsta370set4 (1:14), 'o' ,  k,  errorl0a370set4 (1:14)  ,  ' 

x' ,  k,  errorlsta370set4 (1:14) , k, errorl0a370set4 (1:14) ) 

legend (' xcorr  without  denoising' , ' sta ' ) 

title ('4TH  ORDER  method  with  cross  corellation  after 

denoise;  100  realizations  of  data_370_set4 ' ) 

ylabel ( 'MSE' ) 

xlabel ( 'SNR' ) 

figure ( 7 ) 

plot (1:2047, xy) 

title (' Correlation  Function  of  x  and  y  signals  w/noise') 

figure  (8) 

plot (1:2047, XY) 

title (' Correlation  Function  of  x  and  y  signals  without 
noise ' ) 
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save  errorl0a370set4 ; 
save  errorlsta370set4 

save  Hsta370set4; 
save  H10a370set4; 


MODIFIED  APPROXIMATE  MAXIMUM  LIKELIHOOD  DENOISING 


TECHNIQUE 


<^********************************************************** 

o 

*  *  *  *  *  *  *  *  *  *  *  * 

$.*****************★★************************* *******  ******* 
o 

************ 

%  Denoise:  Approximate  maximum-likelihood  delay  estimation 
%  via  orthogonal  wavelet  transform.  In  this  function, 

%  we  assumed  that  noise  is  Gaussian  noise  and  it  has 


a 

%  flat  freq  response.  We  modify  each  detail  function 

%  by  multipliying  modified  AML  coefficient  based  on 

%  the  signal  and  noise  powers. 

% 

%  SYNTAX:  y=amll (xn, yn) 

o. 

*6 

%  INPUT:  xn  =  Received  signal  from  first  receiver 
%  yn  =  Received  signal  from  second  receiver 

Q. 

O 


%  OUTPUT: 

%  y  =  Denoised  signal 

o 

*o 


%  SUB_FUNC:  None 
%  Written  by  Spiros  Mantis 

o 


function  y=denoise (xn, yn) ; 

xyn=xcorr (xn, yn, ’biased') ; 
[sigmas  b]=max(xyn); 
rx=xcorr(xn,  'biased') ; 
maxx=rx (length (xn) ) ; 
ry=xcorr(yn, 'biased'); 
maxy=ry (length (yn) ) ; 
sigmanl=maxx- sigmas; 
s igman2=maxy- sigmas ; 
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nx=floor (log2 (length (xn) ) ) ; 
ny=floor (log2 (length (yn) ) ) ; 

[cx  lx] =wavedec (xn,nx, 'db4 ' ) ; 

[cy  ly]=wavedec(yn,ny, 'db4 ' ) ; 

dxc= [ ] ; 
for  i=l:nx 

d=detcoef (cx, lx,  i) ; 
dl=length (d) ; 
sigmad= (1/dl) *sum(d.A2) ; 
sigmasd=sigmad-sigmanl ; 
sigmaa (i) = (2A (nx-1) ) *sigmasd; 
if  sigmasd<=0 
wd=0  ; 
else 

wd=sigmasd/ (sigmanl*sigman2+sigmasd* (sigmanl+sigman2 ) ) ; 
end 

dc=wd*d; 
dxc=[dc  dxc] ; 

end 

a=appcoef (cx,lx, 'db4 ' ,nx) ; 
al=length (a) ; 
aux=sum( sigmaa) ; 
sigmasa= ( (2Anx) ) *sigmas-aux; 
if  sigmasa<=0 
wa=0  ; 
else 

wa=sigmasa/ (sigmanl*sigman2+sigmasa* (sigmanl+sigman2 ) ) ; 

end 

ac=wa*a; 
dxc=[ac  dxc] ; 

xd=waverec (dxc, lx, ' db4 ' ) ; 
y=xd; 
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9 .********************************************************** 
o 

'k'k-k'k'k'k’k'k'k'k'k-k 

£.***********************************************★******★*** 

o 

'k'k'k'k-k'k'k-k'k-k'k-k 

% 

%  Modified_tez5 :  This  is  a  test  program  for  modified  AML 
technique.  In  this 

%  program,  we  used  GSM  signal. 

Q, 

"O 

%  SYNTAX:  Modif ied_tez5 
% 

%  INPUT:  None 

o. 

o 

%  OUTPUT:  Mean  square  error  versus  SNR 

o. 

o 

%  SUB_FUN:  Denoise. m 
%  Written  by  Spiros  Mantis 

£.************★********************************************* 

o 

o 

************ 

clear  all 

gsm_set;  %  Configuration  variables  created  in  memory. 
%these  are: 

%  Tb (=  3 . 692e-6) 

%  BT (=  0.3) 

%  OSR (=  4) 

%  SEED (=  931316785) 

%  INIT_L (=  260) 

data=data_gen (INIT_L) ;  %  this  creates  a  binary  data 
[tx_burst, I, Q] =gsm_mod (Tb, OSR, BT, data, TRAINING) ; 

s=I+j*Q; 
sl=length (s)  ; 

pow= (1/sl) *sum (abs (s) . ^2)  ; 


K=100  %  number  of  realizations 

rand ( ' seed' ,40) ; 
f=150*rand(l,K) ; 
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delay=floor (f ) ; 

delay (l:K/2)=-delay (l:K/2)  ;  %  delay  is.  between  -150  to  +150 


n= [20  18  17  16  15  14  13  12  9  6  3  0  -3  -6]; 

SNR=10.A (n./lO) ; 

for  k=l : length (SNR) 

oran=SNR (k) 
for  i=l:K 

x= [zeros (1, 200)  s  zeros (1,224) ] ; 

y= [zeros (1, 200+delay (i) )  s  zeros (1, 224-delay (i) )] ; 
randn ( ' state ' ,  2* (i+j ) ) ; 

noil_real=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 3* (i+j ) ) ; 

noil_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' ,  4* (i+j ) ) ; 

noi2_real=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 5* (i+j ) ) ; 

noi2_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 

noil=noil_real+ j  *noil_imag; 
noi2=noi2_real+j  *noi2_imag; 

xn=x+noil;  %  x  +  noise 
yn=y+noi2;  %  y  +  noise 

%  TDOA  calculation  with  xcorr (  without  denaising) 

xy=xcorr (xn, yn) ;  %  correlation  of  x  and  y  with  xcorr 
% (with  the  presence  of  noise) 

[al  bl] =max (real (xy) ) ; 

erl (i) =delay (i) - (bl-1024 ) ; 

%AML 

x_real=denoise (real (xn) , real (yn) ) ;  %  Denoising  of  the 
%received  signals 
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y_real=denoise (real (yn) , real (xn) ) ; 
x_imag=denoise (imag (xn) , imag (yn) ) ; 
y_imag=denoise (imag (yn) , imag (xn) ) ; 

X=x_real+ j  *x_imag; 

Y=y_real+ j  *y_imag; 

XY=xcorr (X, Y) ;  %  Correlation  of  X  and  Y  (denoised) 

[a2  b2] =max (real (XY) ) ; 

er8 (i) =delay (i) - (b2-1024)  ; 


end 

error8a (k) = (1 /length (er8 ) ) *sum (er8 . A2) ; 
error la (k) = (1 /length (erl) ) * sum (erl . A2) ; 

H8a (k,  : ) =er8 ; 

Hla (k, : ) =erl; 


end 

figure (6) 

k=  [20  18  17  16  15  14  13  12  9  6  3  0  -3  -6] ; 

plot (k, error la (1:14),' o', k,  error 8  a (1:14),' x',k, error la (1:14 
) , k,error8a (1:14)) 

legend (' xcorr  without  denoising ' , ' ami ' ) 

title (' Denoise  AML  method  with  cross  corellation  after 

denoise;  100  realizations') 

ylabel ( 'MSE' ) 

xlabel ( ' SNR' ) 

figure (7) 

plot (1:2047, xy) 

title  (' Correlation  Function  of  x  and  y  signals  w/noise') 

figure (8) 

plot (1:2047, XY) 

title (' Correlation  Function  of  x  and  y  signals  without 
noise ' ) 

save  error8a; 
save  errorla; 

save  Hla; 
save  H8a; 
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TIME-VARYING  MODIFIED  APPROXIMATE  MAXIMUM  LIKELIHOOD 


DENOISING  TECHNIQUE 


%********************************************************** 

************ 

%********************************************************** 

************ 

%  ami 2 :  We  modified  the  MAML  technique  by  dividing  each 
%  detail  function  into  two  segments.  Different 

%  coefficients  for  each  segment  are  computed. 

o. 

'o 

%  SYNTAX:  y=aml2 (xn, yn, delay ) 

% 

%  INPUT:  xn  =  Received  signal  from  first  receiver 

%  yn  =  Received  signal  from  second  receiver 

%  delay  =  True  TDOA  between  xn  and  yn 

%  OUTPUT:  y  =  Error  between  true  TDOA  and  estimated  TDOA 

Q. 

"0 

%  SUB_FUNC:  None 
%  Written  by  Unal  Aktas 

^**************************^r^*^****:*r*^^^:A.^^*1lr^^^^^:A.^^^^^^^^ 

************ 

%********************************************************** 

************ 


function  y=aml2 (xn, yn, delay)  ; 

xyn=xcorr (xn, yn, 'biased' )  ; 
[sigmas  b]=max(xyn); 
rx=xcorr (xn, 'biased' )  ; 
maxx=rx (length (xn) ) ; 
ry=xcorr (yn, 'biased' ) ; 
maxy=ry (length (yn)  )  ; 
sigmanl=maxx-sigmas; 
sigman2=maxy-sigmas; 

nx=floor (log2 (length (xn) )) ;  %10 
ny=f loor (log2 (length (yn) ) ) ; 
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[cx  lx] =wavedec (xn, nx, ' db4 ' ) ; %wavelet  decomposition  of  xn 
%at  level  10 

[cy  ly] =wavedec (yn, ny, ' db4 ' ) ; 

dxc=  [  ]  ; 
for  i=l:nx 

d=detcoef (cx, lx, i) ;  %detail  coefficients  at  level  i 
dl=length (d) ;  %7 

Nsl=128/2A ( i— 1 ) ;  %  the  length  of  the  subblock 

if  Nsl<=l  %i>=7;  l<i>10 

Ns=dl;  %7 
else 

Ns-Nsl;  %  . 25<Nsl>128 

end 

D=ceil (dl/Ns)  ; 
if  dl<Ns*D 

dm=[d  zeros (1, D*Ns-dl) ] ; 

end 

for  k=l:D 

p= ( k-1 ) *Ns+l : k*Ns ; 
sigmad= (1/Ns) *sum(dm(p) . A2) ; 
sigmasd=sigmad-sigmanl ; 
if  sigmasd<=0 
wd=0; 
else 

wd=sigmasd/ (sigmanl*sigman2+sigmasd* (sigmanl+sigman2) ) ; 
and 

dc (p) =wd*dm (p) ; 

end 

dxc= [dc (1 : dl)  dxc] ; 

end 

a=appcoef (cx, lx, ' db4 ' , nx) ; 
al=length (a) ; 
sigmaa= (1/al) *sum(a. A2) ; 
sigmasa=sigmaa-sigmanl; 
if  sigmasa<=0 
wa=0 ; 
else 

wa.=sigmasa/  (sigmanl*sigman2+sigmasa*  (sigmanl+sigman2)  )  ; 

end 

ac=wa*a; 
dxc=[ac  dxc] ; 
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xd=waverec (dxc, lx, 'db4 ' ) ; 


dyc= [ ] ; 
for  i=l:ny 

dy=detcoef (cy, ly, i) ; 
dyl=length (dy) ; 

Nsl=128/2A  (i-1) ;  %  the  length  of  the  subblock 

if  Nsl<=l 
Ns=dyl; 
else 

Ns=Nsl; 

end 

D=ceil (dyl/Ns) ; 
if  dyl<Ns*D 

dmy=[dy  zeros (1, D*Ns) ] ; 

end 

for  k=l:D 

p= (k-1) *Ns+l : k*Ns; 
sigmady= (1/Ns) *sum(dmy (p) .  A2)  ; 
sigmasdy=sigmady-sigmanl; 
if  sigmasdy<=0; 

wdy=0; 

else 

wdy=sigmasdy/ (sigmanl*sigman2+sigmasdy* (sigmanl+sigman2) ) ; 
end 

dcy  (p)  =wdy*dmy  (p) 

end 

dyc= [dcy (1 :dyl)  dye] ; 

end 

ay=appcoef (cy, ly, ' db4 ' , ny) ; 
ayl=length (ay) ; 
sigmaay= (1/ayl) *sum(ay. A2) ; 
sigmasay=sigmaay-sigmanl; 
if  sigmasay<=0 
way=0 ; 
else 

way=sigmasay/ (sigmanl*sigman2+sigmasay* (sigmanl+sigman2 ) ) ; 
end 

acy=way*ay; 
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dyc=[acy  dye]; 
yd=waverec (dye, ly, ' db4 ' ) ; 

% [wdl  wd2  wd3  wd4  wd5  wd6  wd7  wd8  wd9  wc9] 
rxyd=xcorr (xd, yd) ; 

[a5  b5] =max (rxyd) ; 

er5=delay- (b5-1024)  ; 
y=er5; 
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%********************************************++*+*+*****+*+ 

************ 

%****************************+******+*******+***+*+++*+++++ 

************ 

%  tez5t  :  This  is  a  test  program  for  time  varying  MAML 
technique . 

%  In  this  program,  GSM  signal  is  used 

% 

%  SYNTAX:  tez5t 

o 

'o 

%  INPUT:  None 

o 

'o 

%  OUTPUT:  Mean  square  error  versus  SNR 

q. 

"o 

%  SUB_FUNC:  aml2.m 
%  Written  by  Unal  Aktas 

^t*****************************************************^^ 

************ 

%********************************************************** 

************ 


clear  all 


gsm_set;  % 
these  are: 

% 

% 

q. 

■q 


o, 

o 


Configuration  variables  created  in  memory, 

Tb (=  3 . 692e-6) 

BT (=  0.3) 

OSR (=  4) 

SEED (=  931316785) 

INIT  L (=  260) 


data— data_gen ( INIT_L) ;  %  this  creates  a  binary  data 
[tx_burst, I, Q] =gsm_mod (Tb, OSR, BT, data, TRAINING) ; 

%s=I+j*Q; 

data_740_setl 

s=transpose (s2) 

sl=length (s2) ; 

pow= (1/sl) *sum (abs (s) . A2 ) ; 


K=100  %  number  of  realizations 

rand ( ' seed' ,40) ; 
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f=150*rand (1, K) ; 
delay=floor (f ) ; 

delay (1 : K/2) =-delay (1 : K/2) ;  %  delay  is  between  -150  to 

+150 


n= [20  18  17  16  15  14  13  12  9  6  3  0  -3  -6]; 

SNR=1 0 . A ( n . / 1 0 ) ; 

for  k=l : length (SNR) 

oran=SNR (k) 
for  i=l:K 

x= [zeros (1, 259)  s  zeros (1,224) ] ; 

y= [zeros (1, 259+delay(i) )  s  zeros (1, 224-delay (i) )] ; 
randn ( ' state ' , 2* (i+j ) ) ; 

noil_real=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 3* (i+j ) ) ; 

noil_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 4* (i+j ) )  ; 

noi2_real=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 5* ( i+j ) ) ; 

noi2_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 

noil=noil_real+j  *noil_imag; 
noi2=noi2_real+j  *noi2_imag; 

xn=x+noil;  %  x  +  noise 
yn=y+noi2;  %  y  +  noise 


%  TDOA  calculation  with  xcorr (  without  denoising) 

xy=xcorr (xn, yn) ;  %  correlation  of  x  and  y  with  xcorr 
[al  bl] =max (real (xy) ) ; 

erl (i) =delay (i) - (bl-1024)  ; 

%AML 
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el=aml2 (real (xn) , real (yn) , delay (i) ) ; 
e2=aml2 (imag (xn) , imag (yn) , delay (i) ) ; 

ela=amll (real (xn) , real (yn) , delay (i) ) ; 
e2a=amll (imag (xn) , imag (yn) , delay (i) ) ; 

er8 (i) = (el+e2) /2; 

er8a (i) = (ela+e2a) /2; 


end 

error 8 t (k) = (1 /length (er8) ) *sum (er8 . A2) ; 
error la (k) = (1 /length (erl) ) *sum(erl . A2) ; 

error8a (k) = (1 /length (er8a) ) *sum (er8a . A2 ) ; 

H8t (k, : ) =er8 ; 

Hla (k, : ) =erl; 


end 

%load  error8a 
figure (6) 

k= [20  18  17  16  15  14  13  12  9  6  3  0  -3  -6] ; 

plot (k, error la (1 : 14) , ' o' , k, error8t (1:14), ' x' , k, error 8a (1:14 

}  '  d ' 

)  r  u  r  •  •  • 

k, error la (1 : 14) , k, error8t (1 : 14) , k, error8a (1 : 14) ) 

legend ( 'xcorr  without  denoising ' , ' time  varying  ami ' , ' ami ' ) 
title ( 'TVAML  Method  ADS  input  data;  Ts  =  740nsec  Set#l;  100 
realizations ' ) 
figure  (7) 

plot (1 : 541, real (s2) ) 

title ('ADS  signal  (s)  real  component;  740ns') 

save  error8t; 
save  errorla; 

save  Hla; 
save  H8t; 
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APPENDIX  F.  GUIDE  FOR  USING  GSM  SIMULATION  IN  HEWLETT- 


PACKARD  ADVANCED  DESIGN  SYSTEM 


STEP  1:  Start  Program 

If  there  is  no  icon  on  startup  window: 
Start  ->  Programs  =>  HP- ADS  =>  ADS 


STEP  2:  Open  GSM  Example  Project 
From  ADS  (main)  window: 

File  ->  Example  Project  =>  Com_Sys  =>  gsmjprj  =>  ok  =>  wait  a  few  seconds 
for  windows  to  pop  up.  (two  windows  come  up:  GSM_SYS  Schematic^  and  MODEM 
Schematic:  1 .  The  Modem  Schematic  is  for  a  more  detailed  simulation.  I  have  been  primarily 
working  with  the  GSM_SYS  schematic.)  DO  NOT  MAKE  ANY  CHANGES  TO  EITHER  OF 
THESE  FILES! ! ! !  Any  changes  will  be  saved  in  the  system  and  they  will  remain  there  for  future 
users  of  the  program.  Save  to  your  H:  drive  or  other  drive  of  your  choice  and  make  changes  there. 

STEP  3:  Modifying  the  system 

The  options  here  are  unbounded.  You  can  experiment  with  all  the  different  components  in  the 
Component  Library.  But  for  now,  we  just  want  to  know  how  to  get  data  out  of  this  system. 

This  is  done  by  placing  ‘data  sinks’  in  the  desired  locations. 

Go  to  the  common  component  box  at  the  top  of  the  schematic  window  and  click  the 
arrow.  This  will  show  all  the  different  components  available.  => 

For  data  extraction,  you  want  to  use  ‘Sinks.’  From  there  you  can  find  the  appropriate 
‘sink’  for  your  data.  (Timed  and  Spect.  Analyzer  are  the  most  commonly  used)=> 

Click  once  on  the  sink  of  your  choice  => 

Move  the  cursor  into  the  schematic  area  and  place  the  sink  in  the  desired  position  by  left- 
clicking  when  you  have  it  in  a  clear  spot.  => 

End  the  command  by  using  the  right-click  button  on  the  mouse  then  click  ‘end 
command’  or  just  hit  the  escape  key.  => 

Use  wire  to  connect  it  to  the  system,  (there  is  a  wire  button  on  the  third  row  of  toolbar,  or 
you  can  go  to  the  component  button  at  the  top)  => 

Notice  the  name  of  the  sink,  (this  is  listed  under  the  device  type)  This  is  important, 
because  you  will  need  to  recall  this  name  when  you  want  to  plot  or  export  the  data.  You 
can  change  the  name  by  clicking  once  on  it  and  then  typing  the  new  name  right  there  in 
the  same  spot. 
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STEP  4:  Simulations 


After  you  have  your  desired  system  configuration,  are  ready  to  simulate  the  system: 

Simulate(on  top  toolbar)  =>  simulate  =>  wait  for  “simulation  finished”  in  summary  window.  The 
simulation  is  now  complete  and  you  can  extract  your  data. 

STEP  5:  Displaying  Data 

Window(on  top  toolbar)  =>  New  Data  Display  =>  wait  for  data  display  window 
To  see  an  x-y  plot: 

Click  once  on  the  “plot”  icon  =>  move  cursor  into  white  area  of  display(a  box  appears), 
position  the  box  and  left-click.  =>  choose  what  data  you  want  to  plot.  =>  ok  =>  your  plot 
appears  in  the  data  display  window. 

To  see  data: 

Click  once  on  the  “123456”  icon  =>  follow  the  same  procedure  as  above. 

STEP  6:  Extracting  and  Exporting  Data 

I  haven’t  figured  out  how  to  export  the  plots  but  here  is  how  to  extract  the  data: 

In  the  schematic  window:  click  on  “start  the  instrument  server”  button(third  button  third  row  of 
toolbar)  =>  Click  the  following: 

Read/Write:  Write 

Write  to:  File 

File  Format:  Citifile  => 

Use  the  browse  button  to  choose  a  directory  and  file  name.  => 

Choose  desired  data  from  Datasets  box.  => 

Click  ‘write  to  file’  button.  =>  your  data  is  now  saved  in  the  specified  directory. 

You  will  have  to  use  wordpad  or  notepad  to  view  this  data.  Then  you  will  have  to  sort  through  all 
the  data  given  to  get  to  the  specific  data.  This  may  be  annoying,  but  rest  assured  the  data  you  want 
is  in  there  somewhere! ! 


STEP  7:  Importing  your  ADS  data  into  Matlab 

After  you  have  located  your  data  in  wordpad  you  can  cut  and  paste  from  there  into  Matlab. 

If  you  want  to  know  how  many  data  points  you  have  it  is  best  to  use  the  Matlab  Editor/Debugger 
(unless  you  have  a  few  hours  to  spend  counting). 

-  Make  the  data  into  a  vector  in  Matlab  and  use  the  size  command  to  determine  the 
length. 

-  Using  this  method  you  can  cut  it  down  to  the  right  size. 
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